Old Dominion University Desearch Foundation 



/ 




V 


U^'V 

/ 


DEPARTMcNT OF MECHANICAL ENGINEERING AND MECHANICS 
COLLEGE OF ENGINEERING 


OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23508 


/ KJ — ' 0 — 
^ 2 7-5 5 " 


ixT- 

SIMPLE NUMERICAL METHOD FOR 
PREDICTING STEADY COMPRESSIBLE FLOWS 


By 

Ernst von Lavante, Principal Investigator 
and 

N. Duane Nelson, Research Engineer 


Progress Report 

For the period ended November 14, 1986 


Prepared for the 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23665 


Under 

Research Grant NAG-1-633 

Mr. Manuel D. Salas, Technical Monitor 

TAD-Theoretical Aerodynamics Branch 


imSA-Cfi-180345) “"“J Sis“ 

\ Progress Beport, 14 no¥* i^^oo \ 

Ociv.) 42 p Avail: BUS HC G3/02 

{ 


H87-25981 


Unclas 

0063725 


February 1986 


DEPARTMENT OF MECHANICAL ENGINEERING AND MECHANICS 
COLLEGE OF ENGINEERING 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23508 


SIMPLE NUMERICAL METHOD FOR 
PREDICTING STEADY COMPRESSIBLE FLOWS 


By 

Ernst von Lavante, Principal Investigator 
and 

N. Duane Nelson, Research Engineer 


Progress Report 

For the period ended November 14, 1986 


Prepared for the 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23665 


Under 

Research Grant NAG-1-633 

Mr. Manuel D. Salas, Technical Monitor 

TAD-Theoretical Aerodynamics Branch 


Submitted by the 

Old Dominion University Research Foundation 

P.O. Box 6369 

Norfolk, Virginia 23508 


February 1986 


SIMPLE NUMERICAL METHOD FOR 


PREDICTING STEADY COMPRESSIBLE FLOWS 


E. von Lavante 
Old Dominion University 
Norfolk, VA 


and 

N. Duane Mel son** 

NASA Langley Research Center 
Hampton, VA 


Abstract 

A numerical method for solving the isenthalpic form of the governing 
equations for compressible viscous and inviscid flows was developed. The 
method was based on the concept of flux vector splitting in its implicit form. 
The method was tested on several demanding inviscid and viscous configura- 
tions. Two different forms of the implicit operator were investigated. The 
time marching to steady state was accelerated by the implementation of the 
multigrid procedure. Its various forms very effectively increased the rate of 
convergence of the present scheme. High quality steady state results were 
obtained in most of the test cases; these required only short computational 
times due to the relative efficiency of the basic method. 
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I. Introduction 

In recent years, considerable research has been performed towards the 
goal of computing three-dimensional inviscid and viscous transonic flows about 
realistic configurations. The complexity of these flowfields makes their 
numerical prediction very demanding in terms of the capabilities of the 
numerical method and the computer used for the calculation. Numerical methods 
for such problems must achieve high rates of convergence while providing 
results of good quality on reasonably sized computational grids. Computer 
hardware must have sufficient memory to perform the calculations and be fast 
enough to provide reasonable turnaround. 

Many methods for predicting the inviscid compressible flows about 
realistic three-dimensional bodies have been developed. Perhaps the oldest is 
the explicit MacCormack [1] scheme dating back to 1972. Next came the 

implicit (three-factor ADI) method [2,3] using central difference for the 
spatial flux derivatives. The explicit, multistage Runge-Kutta method with 
central differences for the spatial derivatives [4] and multigrid acceleration 
[5] followed and is the method used widely around the world today. 

The newest methods are the Implicit schemes with flux-vector-splitting 
[6-10]. References [6-8] used a full formulation of the Euler equations. 

References [9] and [10] used an isenthalpic formulation which reduces the 
three-dimensional problem to a set of four partial differential equations. 
The energy equation was replaced by an algebraic expression. 

The present effort continues the work by von Lavante and uses the 
isenthalpic assumption in two-dimensions. With the governing equations 
reduced to three partial differential equations, it is necessary to only solve 
3x3 matrices in t ie block tridiagonal system of equations. This requires 
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about half as much work as solving the 4x4 block tri diagonal systems if the 
isenthalpic assumption was not made (9 versus 16 elements). 


Jameson has pointed out the importance of conserving total enthalpy when 
solving the Euler equations (Ref. [113). In his work, care must be taken to 
ensure its conservation. In the present method, with the isenthalpic 
assumption, the conservation of total enthalpy is assured, a priori. 

The isenthalpic assumption is not without its drawbacks. First of all, 
it is limited to steady state calculations since the substantial derivatives 
of the total enthalpy and pressure are related. However, if the pressure is 
'slowly' varying, the isenthalpic equations may be used. Second, for viscous 
calculations, the maximum freestream Mach number is limited to transonic and 
moderate supersonic values due to the requirement of no heat sources or 
sinks. Finally, viscous results can be only considered approximate, since in 
real flows the total enthalpy changes within the boundary layer. 
Notwithstanding these limitations, the present scheme worked well and produced 
good quality results in cases where the flow is steady or slowly varying. 

II. The Equation of Motion 

As noted previously, there is a large class of problems where only steady 
state solutions are of interest. For inviscid flows, the assumptions of 
steady state flow reduces the energy equation to the simple statement that in 
the absence of heat sources and sinks the total enthalpy will remain constant. 
The energy equation is therefore replaced by a simple algebraic equation, 
reducing the number of PDE's to be solved by one. In the case of viscous 
flows, the above statement is not true. However, it is well known that for 
the Prandtl Number Pr=l and adiabatic walls, the total enthalpy will still be 
constant. Many investigators have applied 'e enthalpy damping acceleration 
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technique introduced by Jameson [5] and have effectively driven the total 
enthalpy to zero at the resulting (hopefully) steady state. Several 

investigators used this techique to predict very complex two- and three- 
dimensional configurations and reported results that were in good agreement 
with experimental data. Because of this experience, the present viscous 

formulation assumes that the total enthalpy is constant even in the presence 
of natural dissipation. Due to the obvious limitations of the present 
formulation, only the thin shear layer form of the viscous terms was 

implemented. Here the viscous terms in the normal direction are assumed to be 
much larger than those in the streamwise direction, which is consequently 

neglected. The two-dimensional Navier-Stokes equations for compressible flows 
in vector form for general, body fitted coordinates written in nondimensional 
strong conservation law form using the thin shear layer assumption are 
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Using the definition of the speed of sound at stagnation conditions, 
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resulting in the following form of the equation of state 

p = £ a - (u^ + v^)} (3) 

All variables are nondimensionallzed by the stagnation values (for details, 
see Reference [10]): the primes denoting non-dimensional quantities are 

dropped for convenience. The viscous terms are shown after the application of 
the Stokes hypothesis for bulk viscosity. In the above equations, p Is 

density, u and v are the cartesian velocity components and p Is static 
pressure. The viscous terms Included In Gy will be discussed later. The 
metric coefficients of the transformation of coordinates are defined as 


J X. 


'y- 
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n* 


(4) 


where J Is the Jacobian of the transformation 

J - i/(x^y^ - x^j) (5) 

and U- and U are the contravarlant velocities 
^ h 

U_ = uE + VE U = un + vn. 

E X ^y ti 'x y 

The equivalent Inviscid Euler equations are obtained from Eq. (1) by 
setting Gy ® 0. The development of the solution algorithms for the Euler 
equations follows below. 

III. Development of Inviscid Algorithm 

An Implicit Euler single step temporal scheme was selected for advancing 
the solution of Eq. (1) In time. After linearization In time using the Taylor 
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series expansions of the flux vectors F and G, and approximate factorization 
of the Implicit operator (details are given In Reference [3]), the basic 
algorithm has the form 

[I + At b- a"] [I + At a b"] aq" = -At (a/" + a g") = r" (7) 

gn+1 , qn ^ ^qn 

where A and B are the Jacobian matrices 

The Jacobian matrices A and B are given In detail In Reference [10]. The 
special discretization of Eq. (7) can be carried out In many different ways. 
In the present method, the flux vector splitting approach applied to cell 
centered finite volume formulation was selected. The main reason was Its 
superior ability to capture relatively strong shocks within at most ttio zones. 
It can be also shown that Its truncation error provides the minimum necessary 
damping to limit spurious oscillations In the weak solutions to the Euler 

equations. Based on our previous experience reported In Ref. [10] as well as 

results presented to Ref. [8], It was decided to use the flux vector splitting 
Introduced by van Leer [13] coupled with the so-called MUSCL type 
differencing. The van Leer splitting was selected because the split flux 
vectors are smooth and have smooth first derivatives with respect to the Mach 
number, so that their eigenvalues are also smooth. 

The Inviscid flux vectors F and G each have a complete set of three real 

eigenvectors and can be therefore split Into two vectors, one with non- 

negative eigenvectors and one with non-positive eigenvectors. Following 
Reference [13], these are 
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F = F+ + F~, G = G"^ + G' (9) 

where, for example, F"*" = (f|, F^, F^^, In more detailed, they are in 
cartesian coordinates (denoted here by *) 


• p c.,x 4 (N.,x * ■ P '*.x 4<«*,x - !»' 
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The above split fluxes F and G have, by definition, one zero 

A ^ 

eigenvalue and two positive eigenvalues; similarly, the fluxes F and G 
have one zero and two negative eigenvalues; these are in the case of F: 


'1,2 


(M. 
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\a = 0 

where c* and M* are either and or c*^y and M*^y. 

These split fluxes are, unfortunately, formulated In cartesian 
coordinates only. They have to be transformed into general coordinates 5 
and which is accomplished by simply rotating the local coordinates at a 

given point in the flowfield to a direction parallel with one of the covariant 
vectors t and t[. This procedure is described in some detail in Reference 
[8]; the resulting transformation is 


where 


F » (Q) 




The new dependent variable vector Q is obtained from Q by replacing the 
cartesian velocity components u and v by the physical velocity components u 
and V in the covariant direction t. These are, respectively 


S = (y^ u - x^v)/ + y^ 


; = (x^u ♦ y, V)/ x^ ♦ y^ 
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Knowing the eigenvalues of the split fluxes. It Is now obvious that In 
the special differences In Eq. (7) F'*’ and G'*’ have to be backward differenced 
and F“ and G" have to be forward differenced. This Is accomplished by the 
application of the MUSCL type differencing, described In more detail In 
References [10, 8]. Here, Instead of using the traditional backward or 
forward finite differences operating on F'*’, G"*”, F” and G“, the dependent 
variables Q, which are better differentiable than the flux vectors, are 
extrapolated to the cell faces In positive or negative direction, depending on 
the sign of the eigenvalues. The RHS of Eq. (7) becomes 

r" = - At - •"l-l/2,j ‘" 1 + 1 / 2 , j ■ '" 1 - 1 / 2 , j S,j+l/2 


“ ^1,j-l/2 ^ S,j+l/2 “ S,j-l/2^ 


(14) 


where, for example, 

^1+1/2, j “ ^*^^1+1/2, * *"l+l/2,j “ *" ^^1+1/2, 

'^1,j+l/2 * '^^^'^1, 0+1/2^ ’ S,j+l/2 “ 1,j+l/2^ 

and 

^ 1+1/2, j “ '^1+l,j " '^s 7 ^^1+2, j ‘ ^1+1, 

^1+l/2,j “ ^1,j '^s 7 

etc., with similar expressions In the j-directlon. 

The parameter k^ switches between first order formulation (k^ = 0) and 
second order formulation (kg * 1). 

The present formulation, when applied to transonic and low supersonic 
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flows, did not require the use of flux limiters for essentially oscillation 
free shocks. This was noticed by Anderson, Thomas and van Leer [8] and von 
Lavante and Haertl [10] and was explained in more detail by van Leer [14]. 
The favorable behavour of the present formulation is due to the fact that at 
transonic speeds, the backward running characteristic variable that is being 
extrapolated from downstream of the shock is much smaller than the forward 
running characteristic variable. Despite tte above linear extrapolation, no 
or very small overshoots were encountered. 

The implicit left hand side of Eq. (7) underwent similar modifications as 
the right hand side. The Jacobian matrices A and B were replaced by the 
corresponding Jacobians of the split flux vectors, yielding following form of 
the left hand side of Eq. (7). 

[I + A"^ + fl [I + 5^ + 5^ B“] AQ" = R" (15) 

where A — ■ , A ~ ^ ~ , B ~ gjj and 5^ and 5^ are first 

order backward and forward differences, respectively. The exact form of these 
Jacobian matrices will be given in the full paper. A standard block 
tri diagonal solver was used to solve the system of algebraic equations given 
by Eq. (15). 

An acceleration of the convergence to steady state conditions was 
achieved by the use of local time steps. In this procedure, the time step 
used in each of the cells was determined from the maximum local eigenvalue 
after each iteration. The two-factor, block tridiagonal form of the resulting 
algorithm, given by Eqs. (15) and (14) is relatively easy to vectorize. 

The Eq. (7) is not the only possible Euler implicit form. The 

factorization of the left hand side can be carried out in many different ways. 
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One alternative factorization that seemed to yield more efficient method for 
execution on scalar computers Is a four-factor form 

[I + A"^] [I + A*]CI + sjj B"^] [I + B*] aQ" = r" (16) 

This form of the Implicit operator results In block-bl diagonal L-U 
systems of equations that can be solved very efficiently, each of them In one 
sweep. The Implicit operator based on Eq. (16) was therefore also tested In 
the present Investigation. The directions of sweep were permutated In order 
to preserve the symmetry of the algorithm. However, the L-U scheme was only 
marginally better than the block-tridiagonal scheme, since Its maximum stable 
CFL number was lower, probabily due to Its larger splitting error. Attention 
was therefore focused on the scheme based on Eq. (15), since. In addition to 
being slightly more efficient. It also makes the Inclusion of central 
difference viscous terms much easier than the L-U scheme. 

IV. Development of Viscous Algorithm 

The simplicity of the present viscous scheme Is striking. This Is mainly 
due to the assumption of constant total enthalpy and application of the thin 
shear layer form of the Navler-Stokes equations. Equation (10) already 
Indicated that only two lines have to be added to the flux vector G, In 
addition to the necessary changes In the Implicit operator. 

After some simple manipulations, Gy can be rewritten as 


G.. = |iJ 


-(• 


4 2 


^y- -x^) 

3 *5 U % - 4 >‘5 * \ 


(17) 
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It should be noted that in this case this is exactly the same form as 
obtained by Swanson and Turkel [12] after evaluating the first order 
derivatives in the viscous terms using Green's theorem and then using 
following center difference discretization at the cell face at j+1/2: 

% “ ^,j+l ” 

“ ^i.j ■ ’^i-l.j' 

*“7 ‘“U * 

In the case of the implicit part, the full viscous flux Jacobian matrix 
is relatively complicated and time consuming to evaluate. However, the basic 
underlying assumption of the validity of the present method was that only 
steady state results were required. From the delta form of the governing 
equations, given in Eqs. (17), (15) and (16), it is obvious that if the method 
converges to steady state, the steady state results will not be effected by 
the implicit operator. The dissipative part of the implicit operator can be 
therefore simplified, as long as the stability limits of the scheme are not 
reduced. 

It has been shown previously (References [15, 16] that the implicit part 
of the Navier-Stokes solver (15) with (17) can be significantly simplified by 
substituting the correct Jacobian matrix of the viscous flux by a diagonal 
matrix Iv, where I is the unity matrix and v is the maximum diagonal 
component of the matrix M that is obtained from the thin shear layer Navier- 
Stokes » juations in the following delta form 


% " ^•,j+l ■ ^,0 




(18) 
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i 

[I + a.A] [I + d B + a M d ] aq" = r" (19) 

^ n r) 

V = Ii/(X^ + y^) 

P J T1 

The final form of the present viscous scheme is given by 

[I + + dl A"] [I + a^ + a^ B" + a iv a ] aq" = r” (20) 

^ ^ n T| Tj h 

where, for example, j+1/2 S*i j+1/2 ^ account for 

the viscous effects in the explicit part R*^, The viscous terms are centrally 
j differenced in both the explicit and implicit parts of Eq. (20). 

I 

i 

V. Multi grid 

The Full Approximation Storage (FAS) multigrid scheme (Ref. [1]) must be 
used since the set of equations are nonlinear. A development of the FAS 
scheme is given below. Consider the problem 

ih = f*^ (21) 

where is a nonlinear operator on a grid, G^, with spacing h. The forcing 

function, f is known and is the solution to the problem on the grid with 

spacing h. If we take u^ as an approximation to with an error of 

yh a - u^, 

Eq. (21) can be written as 

Lh (yh + yh) a fh (22) 

u^ is subtracted from both sides of (18) give: 

lh (yh + yh) _ |_h (yh) a fh . i,h (^h) (23) 

If the terms are smooth, the^ can be represented on a coarser grid, G^*^, 
with spacing 2h. The grid 6^*^ ' ’ formed by deleting every other point in 
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6*^. Therefore, G^. Points are eliminated from G^^ to form G^^ and so 

forth to form G®^, G^®^, etc. Written on the coarse grid G^^, Eq, (23) 
becomes 

l2h (j2h Jn ^ y2hj _ |_2h ^2h ^hj , j2h ^^h . ^h ^^hj ^24) 

n n n 

or 

l2h (u2h) , f2h (25) 


where 

f2h , j2h (^h _ lH yhj ^ (j2h ^^hj 

and is the restriction operator. 

Since Eq. (25) is on a coarser grid than Eq. (22), the numerical solution 
for u^^ is much cheaper to obtain because fewer points are involved. Note that 
the operator used on the coarse grid has the same form as the fine grid 
operator, the grid spacing (h and 2h) being the only difference. Once the 
values of u^^ are obtained, the fine grid iterative solution is updated using 
the following equation: 


^“"^New “ ^“"^Old ^2h ‘ ^n^ ^*^"^01 d* 


(26) 


where is the prolongation operator. 


The restriction operator has two forms. One form is used to restrict the 

dependent variables, (u ); i.e. the flow quantities p, pu, and pv. 

For these, the volume weighted average of the values of the function at 

midcells of the four fine grid cells contained in a course grid cell are used 
to set the value of the coarse grid (See Fig. 1). The other form of the 

restriction operator is for the restriction of residuals, CL" (u")]. A 

( 
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simple summation of the residuals over the four fine grid cells composing the 
coarse grid cell Is performed. 

The restriction operations are performed for all Interior points of the 
flow field. At the outer boundaries, only the values of the functions are 
restricted, with no residual restriction. These values are frozen to the fine 
grid values and are not updated on the coarse grids since a lift-correction 
scheme Is used to set the outer-boundary values on the fine grid. The lift- 
correction scheme was found to be less accurate on the coarse grids and tended 
to slow the convergence. At the airfoil surface, the values are not frozen 
and the same boundary condition was used for all the grids. At the wake cut, 
flow values at ghost points were set equal to the flow values from the proper 
points across the wake on all the grids. 

The prolongation operation used In the current work Is a bilinear 
Interpolation In the computation space of the corrections at the four coarse 
grid cells adjacent to the fine grid midcell (see Figure 2). Volume weighting 
Is not used. 

As a first cut, a fixed V-cycle with four grids was used; a fine grid 
with 209x33 cells and three coarse grids with 105x17, 53x9 and 27x5 cells, 
respectively. The program was constructed to allow the number of Iterations 
on each grid between restrictions and prolongation to be controlled by Input. 
Either first-order (kg*0) or second-order (kg=l) approximations can be used. 
In any combination, on each of the grids. Local time stepping was employed on 
each grid with a single reference CFL number controlling all grids. 

Before discussing multigrid results. It Is necessary to define work 
units. Conceptually, a work unit Is the amount of work required to perform 
one fine grid (6^) Iteration, It follows that the work required to perform an 

j 

i. 
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iteration on the grid is 1/4 for two-dimensional problems. For grid 
the work per iteration is 1/16 and for grid G®^ the work is 1/64. To be 
honest, the work required to restrict from one grid to the next must be 
included since a residual calculation is necessary on both the fine and course 
grids. On the conservative side, this can be estimated as the sum of the work 
to perform a fine grid iteration and a coarse grid iteraction, 1.25 for the 
restriction from G*^ to G^^. (The inclusion of the work required for the grid 
transfer is less important if the residual calculation is a small part of the 
update calculation, such as for highly implicit schemes. It is more important 
for explicit schemes.) The work writes for the present multigrid 
computations, expressed as multiples of the single (fine) grid iteration, are 
therefore obtained from 


wu 


"h ^ ^ Zh T "4h TF ^ " 


8h 


1 


where n^...n 3 ^, are the number of iterations on each grid. The work unit count 
was increased in the case of additional iterations during the prolongation 
process. A more precise accounting for the grid transfers would involve 
timing the transfer and then calculating the work required based on the time 
required. In the present study, the conservative method of adding the work 
required to perform fine and a coarse grid iteration is used. This produces 
work which is high by about 15*. 


VI. Results 

The present method was first tested in the single grid version on several 
transonic flow configurations that included internal as well as external 
cases. Only the more interesting test cases will be shown here. 

Inviscid results - The performance of the inviscid (Euler) scheme can be 
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demonstrated on the several standard transonic airfoil test cases. 

In the present work, the NACA 0012 airfoil at various Mach numbers and 
angles of attack was selected. The flow was predicted using a 201x31 C-grid 
that was generated using elliptic grid generation. There were 141 points on 
the airfoil. Following three flow conditions are frequently encountered in 
literature and will be therefore discussed here: 

a) » 0.8, a (angle of attack) = 0 - this supercritical case Is 

particularly well suited to test the ability of a numerical method to perserve 

the symmetry of the flow. The lift coefficient should be zero, while the 

drag coefficient will be nonzero due to the shock. The present method 

predicted the to be 1.41x10"®, a value that Is certainly acceptably close 

to zero. The drag correspond to = 0.0087, which Is In very good resulting 

Mach number and pressure contours are shown In Fig. 3. The shocks were 

captured within at most two zones and are very crisp; the Mach contours are 

very smooth indicating the absence of spurious oscillations. The convergence 

history Is shown In Fig. 4 for the optimum CFL number of 27.5. At this CFL 

number, the spectral radius of the convergence p._ was approximately 0.969, 

sp 

which Is low for a single grid calculation. The residuals shown in Fig. 4 are 
the L-2 norm of the density residual (marked with crosses) and the maximum 
residual of the density. The correct number of supersonic points was reached 
after only 100 Iterations (Fig 4). 

b) M^ = 0,8, a = 1,25° - another supercritical case; It is well suited 
for testing the performance of the boundary conditions, since the lift Is very 
sensitive to their influence. The results obtained from the present method 
were Cj^ = 0.3617 and = 0.0233. These results are in good agreement with 
data published In Reference [17]. The range of best results was given as = 
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0.3632 - 0.3661 and C^j = 0.0229 - 0.0230, achieved on grids that extended up 
to 96 cords from the airfoil!! The comparison with results published by 
Anderson et al. was also favorable; they reported Cj = 0.363 and = 0.0234. 
The corresponding Mach numbe and pressure contours are shown in Fig. 5. The 
shock on the upper surface was again very well captured. The convergence 
history for this case is presented in Fig. 6. Although not as good as in the 
a = 0 case, the spectral radius was still a reasonable p^p = 0.975 at the 
optimum CFL number of 21. The correct number of supersonic points was 
obtained after 189 iterations. 

c) M = 0.63, a = 2° - this case is subcritical. Here, the main 

CO 

difficulty lies, besides the correct lift prediction, in the drag calculation. 
In the absence of shocks, the Cj should be zero. The present scheme computed 
C]^ = 0.3302 and C^j = 0.0006. Both values are in reasonable agreement with 
results reported by Anderson, Thomas and van Leer [8], given as Cj = 0.332 and 
Cjj = 0.0006. The Mach number and pressure contours can be seen in Fig. 7. 
The residual history in Fig. 8 is somewhat surprising. Although the Mach 
number is lower than in the previous cases, the rate of convergence is by far 
the best. This is reflected in the value of the spectral radius 

Pgp = 0.964, obtained at an optimum CFL number of 18. The flow is 
subcritical with no supersonic points. 

Viscous results - The present simplified viscous algorithm was tested on 
several configurations. The consistency of tte method was first investigated 
by grid refinement study done on compressible subsonic and supersonic boundary 
layer flow on a flat plate. The freestream Mach number was 0.5 and the 
reference Reynolds number was Reg = 5,000. The computations were carried out 
on three grids: 51x51, 51x76 and 51x101 that were exponentially stretched in 

the direction normal to the plate. The grid was refined in the direction 
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normal to the plate. The results will be presented in the full version of the 
paper; they can be summarized by stating that the skin friction coefficients 
as well as the velocity profiles consistently improved as the grid was refined 
and agreed closely with the Blasius solution. 

The viscous method was subseqently tested on a supersonic diffusor 
configuration with inflow Mach number M*2.0 and upper wall compression corner 
that generated on oblique shock. The shock was reflected off the lower wall 
at a Reynolds number Re^^ = 3 10® based on the length along the lower wall, and 
was strong enough to cause separation of the boundary layer. The configura- 
tion is similar to tht used by Thomas and Walters [20]. Computations using 
the present method were done on two grids; 51x51 (coarse) and 51x10 (finer). 
The resulting ratio of static pressure to total pressure as well as skin 

friction coefficient Cf at the lower are shown in Fig. 9 as compared with 

experimental data given by Hakkinen et al. [19]. It can be seen that the 

coarse grid had insufficient resolution to correctly predict the extend of the 

separated zone and the formation of the pressure plateau typical for lambda 
shock structure. Refining the grid resulted in much better prediction of the 
pressure and C^; except for the reattachment region, the agreement with 
experimental results is quite good. This tendency is consistent with results 
given in Ref. [20]; minor differences in the pressure profile can be explained 
by the presence of the upper wall in present configuration. The rate of 
convergence, given in Fig. 10, was rather high. At an maximum local CFL 
number of 1450 (global time steps were used in this case) the maximum residual 
was reduced eight orders of magnitude within 500 iterations. 

More recently, results were obtained for laminar flows over transonic 
airfoil NACA 0012, showir separated flow near the trailing edge. These will 
be included in the full vi 'sion of this paper. 
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Mul tigrid - All the multigrid computations were carried out on a 32 bit 

computer, as compared to the single grid cases shows above, which were 

executed on a 60 bit machine. To show the acceleration of the multigrid 

scheme, results for the test case of a NACA0012 at M = . 80, a * 1.25 are 

<0 

presented. In Figure 11, is the convergence history of a baseline, no 

mul tigrid calculation at optimum CFL number. This calculation was performed 
on a 201x31 cell grid, somewhat coarser than the 209x33 cell grid used for the 
mul tigrid calculations. One hundred and sixty-five iterations were required 
to converge the lift coefficient to within 1 % of the converged solution and 
over 250 iterations were required to reduce the maximum residual three orders 
of magnitude, (see case 1, Table 2). Figure 3 shows the logarithm of the 
maximum residual divided by the initial residual, the lift coefficient, and 
the number of supersonic points are plotted versus number of iterations. 

The baseline calculation was then accelerated using mul tigrid with four 
grids in a V-cycle. As a first cut, two iterations were performed on each of 

the grids for 40 cycles. The convergence history of this case (case 2 of 

Table 10 is shown in Figure 12. The lift converged (within 1 %) in 107 work 
units. In case 3, more iterations were performed on the coarser grids and the 
lift was obtained in 88 work units. The residual also converged more quickly 
(Figure 13). 

Other researchers have successfully used mul tigrid without performing 
iteractions between the prolongation operations (Ref. 5, 8). This was tried 
in case 4 (Figure 14) and was found to reduce the time needed to achieve the 
lift coefficient (83 work units) but produced a slower reduction of the 
maximum residual. 

Since the computer program used in the presen work allows a choice 
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between first and second-order differencing, a scheme with second-order 
differencing on the fine grid and first-order differencing on the coarse grids 
was tried. It was hoped the first-order differencing would provide better 
smoothing on the coarse grids. This was not the case. (See Figure 15). The 
lift required 105 work units to achieve a level within 1 % of the converged 
answer. 

It should be noted that same Investigation (Ref. 17) define work units 

simply as cycles. Using this definition of convergence and acceleration, the 

present method converged in 9-30 cycles, starting from uniform flow. This 
corresponds to an acceleration ratio of 10-12. 

To date, the optimum V-cycle FAS multigrid strategy found with this 
method of solving the 2-D isenthalpic Euler equations is to use second-order 
differencing or all the grids and to perform more iterations on the coarser 
grids than on the fine grid. To most quickly obtain the lift, no iterations 
should be performed between prolongations. It was found that the optimum CFL 
number with the multigrid acceleration was close to the optimum for the single 
grid calculation. 

The present scheme was also used to predict flows about other AGARD 

airfoil configurations at following condition: » 0.85, a = 1°; 

= 1.2, a = 0° and = 1.2, a = 1®. The full version of this paper will 
also include a study of inviscid! separation on a backward facing step, carried 
out with the present scheme. 

The present algorithm is oirrently being extended to three dimensions. 
Some three-dimensional results witri be also Included in this paper. 

i 
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Future Work 

The acceleration obtained by increasing the number of iterations on the 
coarse grids suggested that a W-cycle may yield benefits over the V-cycle 
currently used. This will be investigated. To date, a simple reference CFL 
number from grid to grid should be explored. 
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Table 1 - Summary of Airfoil Results 
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Table 2 - Summary of Results, Multigrid 
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- Schematic of fine acad coarse grid cells used for restriction. 
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Figure 2.- Schematic of fine and coarse grid cells used for prolongation 















Mach nunber contours, pressure contours and sonic lines for NACA 0012 
airfoil at M = 0.8, angle of attack = 1.25°. 








b) Finer grid (51x101) 


■ experirasnt 
^ Blasius solution 
— present results 


9 — Results for viscous diffusor calculations. Skin friction coefficient 
and pressure ratio as oontpared with ej^^erimental da t a. 
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Figure 12 


Pour grid nniltigrid iteration-history (60 cycles) at CFL 
for " 0.8 and a =* 1.25® (Case 2)- 
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Figure 13.- Pour grid multigrid 
for Me, « 0.8 and a 
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Figure 15.- Four grid multigrid iteration-history (30 cycles) at 
. for = 0.8 and a « 1.25° (Case 5). 
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